l# ==================== 加载包 ====================
library(gtsummary)
library(dplyr)
library(flextable)

# ==================== 1. 定义变量列表 ====================
categorical_vars <- c(
  "性别", "吸烟", "饮酒", "肿瘤类型", "肿瘤部位", 
  "T", "N", "辅助治疗", "手术类型", "营养管", 
  "合并症", "并发症", "CNP", "营养不良风险"
)

continuous_vars <- c(
  "体重", "年龄", "三个月体重下降", "住院时间", "出血量",
  "血红蛋白", "红细胞", "血小板",
  "中性粒", "淋巴",
  "白细胞", "白蛋白",
  "总蛋白", "尿素氮", "肌酐",
  "PNI", "PLR", "NLR", "SII", "BMI", "NRS2002", "GNRI"
)

group_col <- "分类"   # 分组变量名，取值为 1 和 3

# ==================== 2. 数据检查与清洗 ====================
if (!exists("qws")) {
  stop("错误：数据框 'qws' 不存在，请先读入数据并赋值给 qws。")
}
if (!is.data.frame(qws) || nrow(qws) == 0) {
  stop("错误：数据框 'qws' 为空或不是数据框。")
}
if (!group_col %in% names(qws)) {
  stop("错误：数据中不存在名为 '", group_col, "' 的列。")
}

cat("分组列 '", group_col, "' 的原始取值分布：\n")
print(table(qws[[group_col]], useNA = "ifany"))

# 过滤分组缺失值
qws_filtered <- qws %>%
  filter(!is.na(!!sym(group_col)))
if (nrow(qws_filtered) == 0) {
  stop("错误：分组列全部为缺失值，无法生成表格。")
}

# 仅保留分组取值为 1 和 3 的行（若存在其他值则过滤）
valid_levels <- c(1, 3)
if (!all(qws_filtered[[group_col]] %in% valid_levels)) {
  warning("分组列含有除 1 和 3 之外的值，将只保留 1 和 3。")
  qws_filtered <- qws_filtered %>%
    filter(!!sym(group_col) %in% valid_levels)
  if (nrow(qws_filtered) == 0) {
    stop("错误：过滤后无有效分组（1 或 3），请检查数据。")
  }
}
qws <- qws_filtered

# 检查变量名是否都存在，若缺失则给出警告并从列表中移除
missing_vars <- setdiff(c(categorical_vars, continuous_vars), names(qws))
if (length(missing_vars) > 0) {
  warning("以下变量在数据中不存在，将被忽略：", 
          paste(missing_vars, collapse = ", "))
  categorical_vars <- intersect(categorical_vars, names(qws))
  continuous_vars <- intersect(continuous_vars, names(qws))
}

# ==================== 3. 数据转换 ====================
# 分组转因子（标签可自定义）
qws[[group_col]] <- factor(qws[[group_col]], 
                            levels = c(1, 3), 
                            labels = c("数据集1", "数据集3"))

# 分类变量转因子
if (length(categorical_vars) > 0) {
  qws <- qws %>%
    mutate(across(all_of(categorical_vars), as.factor))
}

# 连续变量转数值（安全转换，避免因子误转）
if (length(continuous_vars) > 0) {
  qws <- qws %>%
    mutate(across(all_of(continuous_vars), ~ as.numeric(as.character(.))))
}

# 检查分组后是否至少有两个水平
if (length(unique(qws[[group_col]])) < 2) {
  stop("错误：分组列只有一个水平，无法生成分列表。")
}

# ==================== 4. Shapiro-Wilk 正态性检验 ====================
# 存储结果
normality_results <- list()   # 每个变量的 p 值
statistic_list <- list()      # 每个连续变量对应的统计量字符串

if (length(continuous_vars) > 0) {
  for (var in continuous_vars) {
    vec <- qws[[var]]
    vec_clean <- vec[!is.na(vec)]
    n <- length(vec_clean)
    
    if (n >= 3 && n <= 5000) {
      p_val <- shapiro.test(vec_clean)$p.value
      normality_results[[var]] <- p_val
      if (p_val >= 0.05) {
        statistic_list[[var]] <- "{mean} ({sd})"
        cat(sprintf("变量 '%s' 符合正态分布 (p = %.4f) → 使用均数±标准差\n", var, p_val))
      } else {
        statistic_list[[var]] <- "{median} [{p25}, {p75}]"
        cat(sprintf("变量 '%s' 不符合正态分布 (p = %.4f) → 使用中位数[P25,P75]\n", var, p_val))
      }
    } else {
      warning(sprintf("变量 '%s' 有效样本量 %d 不在 Shapiro-Wilk 检验允许范围 (3~5000)，默认使用中位数", var, n))
      statistic_list[[var]] <- "{median} [{p25}, {p75}]"
      normality_results[[var]] <- NA
    }
  }
}

# 打印正态性检验汇总表（便于核对）
cat("\n========== Shapiro-Wilk 检验结果汇总 ==========\n")
if (length(normality_results) > 0) {
  df_norm <- data.frame(
    变量 = names(normality_results),
    p值 = unlist(normality_results),
    正态性 = ifelse(unlist(normality_results) >= 0.05 & !is.na(unlist(normality_results)), "是", 
                  ifelse(is.na(unlist(normality_results)), "未检验", "否"))
  )
  print(df_norm)
} else {
  cat("无连续变量可供检验。\n")
}
cat("==================================================\n\n")

# ==================== 5. 构建统计量列表（仅连续变量） ====================
stat_list <- list()
for (var in continuous_vars) {
  if (!is.null(statistic_list[[var]])) {
    stat_list[[var]] <- statistic_list[[var]]
  }
}

# ==================== 6. 构建两组对比的检验方案 ====================
test_list <- list()

# 分类变量：使用卡方检验（若期望频数<5，可改为 fisher.test）
for (var in categorical_vars) {
  if (var %in% names(qws)) {
    test_list[[var]] <- "chisq.test"
  }
}

# 连续变量：根据正态性结果选择 t.test 或 wilcox.test
for (var in continuous_vars) {
  if (var %in% names(qws)) {
    if (!is.na(normality_results[[var]]) && normality_results[[var]] >= 0.05) {
      test_list[[var]] <- "t.test"
    } else {
      test_list[[var]] <- "wilcox.test"
    }
  }
}

# ==================== 7. 生成含 P 值的基线表 ====================
table1 <- qws %>%
  tbl_summary(
    by = group_col,
    statistic = stat_list,                    # 仅连续变量，分类使用默认
    type = list(all_of(continuous_vars) ~ "continuous"),   # 强制类型，避免误判
    digits = list(
      all_continuous() ~ c(1, 1),            # 连续变量保留1位小数
      all_categorical() ~ c(0, 1)            # 分类变量百分比保留1位小数
    ),
    missing = "no"                            # 不显示缺失值计数
  ) %>%
  add_n() %>%                                 # 添加每行有效例数
  add_p(
    test = test_list,
    pvalue_fun = ~ style_pvalue(.x, digits = 3)   # P 值保留3位小数，<0.001自动显示
  )

# 预览表格（此时已包含分组统计、N、P值）
print(table1)

# ==================== 8. 导出为 Word 文档 ====================
ft <- as_flex_table(table1)

# 自定义列标题（中文）
ft <- set_header_labels(ft,
                         stat_1 = "数据集1",
                         stat_2 = "数据集3",
                         label = "变量",
                         p.value = "P值")

# 保存
save_as_docx(ft, path = "Table1_baseline_with_p.docx")

cat("\n表格已成功导出至：", getwd(), "/Table1_baseline_with_p.docx\n")